/*
 *  distance.c
 *  Sandbox
 *
 *  Created by Joshua Lewis on 10/28/10.
 *  Parallel algorithms for calculating distance matrices.
 *
 */

#include "distance.h"

 // Calculate distances with vector/vector operations (BLAS1)
void distance1(int N, int D, float *data, float *result) {
	float norm;
	int i, j, m, o;
	
	if(1) { // Toggle parallel for testing
		int blockSize = 20;
#pragma omp parallel private(i, j, o, norm)
		{
			float sum;
			float *copy = (float *)malloc(sizeof(float) * D); // Unstable
#pragma omp for schedule(guided)
			for(m = 0; m < N / blockSize; m++)
				for(o = m; o < N / blockSize; o++)
					for(i = m * blockSize; i < (m + 1) * blockSize - (m == o ? 1 : 0); i++)
						for(j = (m == o ? i + 1 : o * blockSize); j < (o + 1) * blockSize; j++)
							if(1) // Toggle BLAS for testing
							{
								cblas_scopy(D, &data[j * D], 1, copy, 1);				
								cblas_saxpy(D, -1, &data[i * D], 1, copy, 1);
								norm = cblas_snrm2(D, copy, 1);
								result[utndidx(i, j)] = norm;
							} else {
								sum = 0.f;
								for(int z = 0; z < D; z++)
									sum += (data[i * D + z] - data[j * D + z]) * (data[i * D + z] - data[j * D + z]);
								result[utndidx(i, j)] = sum;
							}
			free(copy);
		}
	}
	else {
		float *copy = (float *)malloc(sizeof(float) * D);
		
		for(int i = 0; i < N - 1; i++)
			for(int j = i + 1; j < N; j++)
			{
				cblas_scopy(D, &data[j * D], 1, copy, 1);				
				cblas_saxpy(D, -1, &data[i * D], 1, copy, 1);
				norm = cblas_snrm2(D, copy, 1);
				result[utndidx(i, j)] = norm;
				//printf("%f ", norm);
			}
		
		free(copy);
	}
	
	return;
}

// Calculate distances with matrix/matrix operations (BLAS3)
void distance3(int N, int D, float *data, float *result) {
	int i, j, m, o;
	int blockSize = 256;
  int threadNum = omp_get_max_threads(); // Number of threads OpenMP will spawn

	float *diag = (float *)malloc(sizeof(float) * N);
  float *C = (float *)malloc(threadNum * blockSize * blockSize * sizeof(float));
  
#pragma omp parallel private(i, j, o)
	{ 
    int th_id = omp_get_thread_num();
    
#pragma omp for schedule(guided)
		for(m = 0; m < N / blockSize; m++)
		{
			cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, blockSize, blockSize, D,
                  1, &data[m * blockSize * D], D,
                  &data[m * blockSize * D], D, 0,
                  &C[th_id * blockSize * blockSize], blockSize);
			
			for(i = 0; i < blockSize; i++)
				diag[m * blockSize + i] = C[th_id * blockSize * blockSize + i * (blockSize + 1)];
      
			for(i = 0; i < blockSize; i++)
				for(j = i + 1; j < blockSize; j++)
					result[utndidx(i + m * blockSize, j + m * blockSize)] = \
            sqrt(diag[i + m * blockSize] + diag[j + m * blockSize] - \
            2 * C[th_id * blockSize * blockSize + j * blockSize + i]);
    }
		
#pragma omp for schedule(guided)
		for(m = 0; m < N / blockSize; m++)
			for(o = m + 1; o < N / blockSize; o++)
			{
				cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, blockSize, blockSize, D,
                    1, &data[m * blockSize * D], D,
                    &data[o * blockSize * D], D, 0,
                    &C[th_id * blockSize * blockSize], blockSize);
				
				for(j = 0; j < blockSize; j++)
					for(i = 0; i < blockSize; i++)
						result[utndidx(j + m * blockSize, i + o * blockSize)] = \
              sqrt(diag[i + o * blockSize] + diag[j + m * blockSize] - \
              2 * C[th_id * blockSize * blockSize + i * blockSize + j]);
			}
  }
  
  free(C);
	free(diag);

	return;
}